function [vol, t_eq, time, C_tot, C, C_dam, A, Ta_tum, T_lym, Ta_lym, Tb_lym] = radioimmuno_response_model(param, delta_t, free, t_f1, t_f2, D, t_rad, t_treat_c4, t_treat_p1, LQL, activate_vd, use_Markov)

%INFO: numerical code to solve the radioimmunotherapy response model

%% Notation

    % There are two physical locations: tumor and lymphatic activation site

    % TUMOR populations:
    % C: cancer cells
    % C_dam: RT damaged cancer cells "doomed" (Cd in the article)
    % Ta_tum: activated T-cells (Ta in the article)

    % Lymphatic ACTIVATION SITE
    % A: antigens/APCs (hat{A} in the article)
    % T_lym: available T-cells to be activated (hat{T} in the article)
    % Ta_lym: active T-cells (hat{Ta} in the article)
    % Tb_lym: blocked/inactivated T-cells (hat{Tb} in the article)

%% OUT:

    % vol: (int) tumor volume (tumor cells & T-cells)
    % t_eq: (real) day in which initial tumor volume or initial time is achieved
    % time: (vec) temporal discretization vector
    % C_tot: (vec) total tumor cells
    % C: (vec) non-damaged tumor cells
    % C_dam: (vec) damaged tumor cells
    % A: (vec) antigens in T-cell activation area
    % Ta_tum: (vec) activated T-cells in tumor area
    % T_lym: (vec) T-cells available to be activated in the activation zone
    % Ta_lym: (vec) activated T-cells
    % Tb_lym: (vec) inactivated T-cells

%% IN:

    % param: (vec) parameters
    % delta_t: (real) time step
    % free: (vec) controls free growth and start of treatment condition
    %    free(1): (int) 1 for free behavior until an specific time or volume, 0 otherwise 
    %    free(2): (int) 1 for time, 0 for volume
    %    free(3): (real) initial volume or time
    % t_f1: (real) maximum time allowed to reach the initial state
    % t_f2: (real) simulation time after initial state is reached
    % D: (vec) RT dose
    % t_rad: (vec) RT days
    % t_treat_c4: (vec) IT days with anti-CTLA-4
    % t_treat_p1: (vec) IT days with anti-PD-(L)1
    % LQL: =1 if LQL, otherwise (modified) LQ (with param(33)) 
    % activate_vd: =1 if vascular death (D > 15)
    % use_Markov: =1 if Markov, otherwise deterministic

%% param:

    %  1  C_0: (real) initial number of tumor cells, fixed to 100
    %  2  lambda: (real) exponential tumor growth lambda*C*(1-lambda_2*C)
    % 33  lambda_2: (real) modulates tumor growth
    %  3  alpha_C: (real) tumor cells death LQ linear parameter
    %  4  beta_C: (real) tumor cells death LQ quadratic parameter
    % 34  beta_2: (real) tumor cells death LQ quadratic parameter related to the RT dose
    %           In the article it is named c or x depending on the chosen damage model
    %  5  phi: (real) tumor cells RT death rate
    %  6  tau_dead_1: (real) tumor cells death delay, fixed to 1
    %  7  tau_dead_2: (real) tumor cells death delay, fixed to 1.5
    %  8  empty, not used
    %  9  vol_C: (real) tumor cell volume
    % 10  A_0: (real) initial number of antigens
    % 11  rho: (real) antigens natural liberation rate
    % 12  psi: (real) antigens RT liberation rate
    % 13  sigma: (real) antigens natural 'death'
    % 14  tau_1: (real) antigens delay
    % 15  Ta_tum_0: (real) initial number of T-cells in the tumor area
    % 16  alpha_T: (real) T-cell death LQ linear parameter, fixed to 0.1
    % 17  beta_T: (real) T-cell death LQ quadratic parameter, fixed to 0.01
    % 18  tau_2: (real) T-cell delay
    % 19  eta: (real) T-cell natural 'death'
    % 20  T_lym_0: (real) initial number of available T-cells in the activation zone
    % 21  empty, not used
    % 22  h: (real) T-cell proliferation rate
    % 23  iota: (real) T-cell immune elimination rate
    % 24  vol_T: (real) T-cell volume
    % 25  c4: (real) anti-CTLA4 concentration
    % 26  r: (real) activated/non-activated T-cells relational term, fixed to 5
    % 27  ni: (real) anti-CTLA4 consumption rate, fixed to 0.1
    % 28  a: (real) activated T-cell rate
    % 29  p: (real) immune tumor cell death rate
    % 30  q: (real) immune death related to effector/target rate
    % 31  s: (real) steepness of the curve
    % 32  recovery: (real) recovery of vascular damage, fixed to 0.05 
    % 35  p1: (real) anti-PD(L)1 concentration
    % 36  mi: (real) anti-PD(L)1 consumption rate, fixed to 0.5
    
C_0 = param(1);
lambda = param(2);
lambda_2 = param(33);
alpha_C = param(3);
beta_C = param(4);
beta_2 = param(34);

phi = param(5);
tau_dead_1 = param(6);
tau_dead_2 = param(7);
vol_C = param(9);

A_0 = param(10);
rho = param(11);
psi = param(12);
sigma = param(13);
tau_1 = param(14);

Ta_tum_0 = param(15);
alpha_T = param(16);
beta_T = param(16)/10;%param(17);
tau_2 = param(18);
eta = param(19);
T_lym_0 = param(20);
h = param(22);
iota = param(23);
vol_T = param(24);
c4 = param(25);
r = param(26);
ni = param(27);
a = param(28);
b = (r - 1) * a;

p = param(29);
q = param(30);
s = param(31);
p1 = param(35);
mi = param(36);

recovery = param(32);

% Temporal discretization
time = 0: delta_t: t_f1 + t_f2 + 1;
m = length(time);

%% Selection LQL or modified LQ
if LQL == 1 && D(1) > 0
    beta_C = min(beta_C, 2 * beta_C * (beta_2 * D(1) - 1 + exp( -beta_2 * D(1))) / (beta_2 * D(1))^2);
else
    beta_C = beta_C * (1 - beta_2 * sqrt(D(1)));
end

%% Vascular death
if activate_vd == 1 && D(1) > 15
    vascular_death = 0;
else
    vascular_death = 1;
end

%% Variable initialization
C = zeros(1, m);       % Tumor cells (tumor)
A = zeros(1, m);       % Antigens (activation zone)
Ta_tum = zeros(1, m);  % Activated T-cells (tumor)
T_lym = zeros(1, m);   % T-cell available to be activated (activation zone)
Ta_lym = zeros(1, m);  % Activated T-cells (activation zone)
Tb_lym = zeros(1, m);  % Inactivated T-cells (activation zone)
vol = zeros(1, m);     % Tumor volume

% Delay index
del_1 = max(1, round(tau_1/delta_t));
del_2 = max(1, round(tau_2/delta_t));

% Auxiliar variable to take account of tumor cell killing
d = length(D);
C_dead = zeros(1, d);  % Damaged tumor cells at each RT time
M = zeros(d, m);       % Alive damaged tumor cells evolution for each RT dose
C_dam = zeros(1, m);   % Total alive damaged tumor cells at each time step
C_tot = C;             % Total alive tumor cells

% Surviving fraction with LQ model parameters
SF_C = zeros(1, d);    % Tumor cells surviving fraction
SF_T = zeros(1, d);    % Lymphocytes surviving fraction

% Variables initial value
C(1) = C_0;
A(1) = A_0;
Ta_tum(1) = Ta_tum_0;
T_lym(1) = T_lym_0;
C_tot(1) = C_0;

% Free behavior in time or volume
free_flag = free(1);   % 1 for free behavior, 0 otherwise
free_op = free(2);     % 1 for time, 0 for volume

t_eq = -1;
vol_flag = 1;          % 1 if initial volume was achieved, 0 otherwise
time_flag = 1;         % 1 if initial time was achieved, 0 otherwise

if free_flag == 1
    if free_op == 0
        vol_in = free(3);
        vol_flag = 0;
    else
        t_in = free(3);
        time_flag = 0;
    end
else
    m = t_f2/delta_t + 1;
end

p_1 = 0;
c_4 = 0;

tf_id = max(1, round(t_f2 / delta_t));

k = 1;                 % Radiation vector index
ind_c4 = 1;            % c4 treatment vector index
ind_p1 = 1;            % p1 treatment vector index

for i = 1:max(del_1, del_2) + 1
    C(i) = C_0;
    A(i) = A_0;
    Ta_tum(i) = Ta_tum_0;
    T_lym(i) = T_lym_0;
    Ta_lym(i) = 0;
    Tb_lym(i) = 0;
    C_tot(i) = C_0;
end

%% Algorithm

j = i;
im_death = Ta_lym;

while j <= m-1
    
    % Tumor cell
    prol = growth(lambda, C_tot(j), lambda_2);

    p_11 = p_1;
    ind_p11 = ind_p1;

    [im_death(j), p1_flag, p_1] =immune_death_dePillis(C_tot(j), Ta_tum(j), p, q, s, p1, p_1, mi, vol_flag, time_flag, time(j+1), t_treat_p1(ind_p1), delta_t);

    if p1_flag == 1
        ind_p1 = min(ind_p1 + 1, length(t_treat_p1));
    end
    
    %Markov
    if C(j) <= 1000 && use_Markov
        C(j+1) = markov_TCP_analysis(im_death(j), prol, C(j), delta_t);
    elseif C(j) == 0
        %Do nothing
    else
        C(j+1) = max(0, C(j) + delta_t * (prol - im_death(j)) * C(j));
    end

    % T-cells
    [T_lym(j+1), A(j+1), Ta_lym(j+1), Tb_lym(j+1), c4_flag, c_4] = ...
    A_activate_T(a, b, T_lym_0, h, c4, c_4, ni, t_treat_c4(ind_c4), time(j+1), delta_t, T_lym(j), A(j), vol_flag, time_flag);
    
    if c4_flag == 1
        ind_c4 = min(ind_c4 + 1, length(t_treat_c4));
    end

    % Antigens
    nat_rel = natural_release(rho, C_tot( (j+1) - del_1 ));

    %Dead tumor cells releasing antigens
    dead_step = M(:, j-del_1) - M(:, j+1-del_1);
    dead_step(dead_step < 0) = 0;
    dead_step = sum(dead_step);

    RT_rel = RT_release(psi, dead_step);
    A_nat_out = A_natural_out(sigma, A(j));
    A(j+1) = A(j+1) + delta_t * (nat_rel + A_nat_out) + RT_rel;
    
    % T-cell
    T_out = immune_death_T(iota, C(j) + C_dam(j), Ta_tum(j));
    T_nat_out = T_natural_out(eta, Ta_tum(j));
    
    %Vascular death
    Ta_tum(j+1) = Ta_tum(j) + vascular_death * Ta_lym((j+1) - del_2) + delta_t * (T_out + T_nat_out);
    
    if (time(j+1) > t_eq && activate_vd == 1 && D(1) >= 15)
        vascular_death = min(1, recovery * (time(j+1) - t_eq));
    end
    
    if vol_flag == 1 && time_flag == 1 && D(1) ~= 0 && abs(time(j+1) - t_rad(k)) <= delta_t/2
        SF_C(k) = exp(- alpha_C * D(k) - beta_C * D(k) ^ 2);
        SF_T(k) = exp(- alpha_T * D(k) - beta_T * D(k) ^ 2);
        
        C_dead(k) = (1 - SF_C(k)) * C(j+1);
        C(j+1) = C(j+1) - C_dead(k);
        
        Ta_tum(j+1) = Ta_tum(j+1) - (1 - SF_T(k)) * Ta_tum(j+1);
        
        for ii = 1:d
            C_kin = tum_kinetic(phi, tau_dead_1, tau_dead_2, time(j+1) - t_rad(ii));
            im_death_d = immune_death_dePillis(C_tot(j), Ta_tum(j), p, q, s, p1, p_11, mi, vol_flag, time_flag, time(j+1), t_treat_p1(ind_p11), delta_t);            
            M(ii,j+1) = max(0, M(ii,j) + delta_t * (C_kin - im_death_d) * M(ii,j));
        end
        
        M(k,j+1) = M(k,j+1) + C_dead(k);
         
        % The sum of the columns of M is the total damaged tumor cells that
        % are going to die in each time step
        C_dam(j+1) = sum(M(:, j+1));
        
        k = min(k + 1 , length(t_rad));
       
    elseif vol_flag == 1 && time_flag == 1 && D(1) ~= 0
        for ii = 1:d
            C_kin = tum_kinetic(phi, tau_dead_1, tau_dead_2, time(j+1) - t_rad(ii));
            im_death_d = immune_death_dePillis(C_tot(j), Ta_tum(j), p, q, s, p1, p_11, mi, vol_flag, time_flag, time(j+1), t_treat_p1(ind_p11), delta_t);                       
            M(ii,j+1) = max(0, M(ii,j) + delta_t * (C_kin - im_death_d ) * M(ii, j));
        end
        C_dam(j+1) = sum(M(:,j+1));
    end
    
    % Non-negative values control
    if C(j+1) < 0; C(j+1) = 0; end
    if C_dam(j+1) < 0; C_dam(j+1) = 0; end
    if A(j+1) < 0; A(j+1) = 0; end
    if Ta_tum(j+1) < 0; Ta_tum(j+1) = 0; end
    
    C_tot(j+1) = C(j+1) + C_dam(j+1);

    vol(j+1) = tumor_volume( C_tot(j+1), Ta_tum(j+1), vol_C, vol_T );
    
    if vol_flag ~= 1 && vol(j+1) >= vol_in
        t_eq = time(j+1);
        t_rad = t_rad + t_eq;
        t_treat_p1 = t_treat_p1 + t_eq;
        t_treat_c4 = t_treat_c4 + t_eq;
        vol_flag = 1;
    elseif time_flag ~= 1 && time(j+1) >= t_in
        m = j + 1 + tf_id;
        t_rad = t_rad + t_in;
        t_treat_p1 = t_treat_p1 + t_in;
        t_treat_c4 = t_treat_c4 + t_in;
        time_flag = 1;
    else
        j = j + 1;
    end   

    if time(j-1) > t_f1 && vol_flag == 0, return, end
    if time(j-1) > t_eq + t_f2 && vol_flag == 1, return, end
end %while

end %function


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                      %
%      Tumor cells natural growth      %
%                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = growth(lambda, C_tot, lambda_2)
    % Logistic proliferation
    f = lambda * (1 - lambda_2 * C_tot);
end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                      %
%      Immune death controlled by PD-L1/PD-1 axis      %
%                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [f, m, p_1] = immune_death_dePillis(C, T, p, q, s, p1, p_1, mi, vol_flag, time_flag, t, t_treat, delta_t)
    m = 0;
    if vol_flag == 0 || time_flag == 0
    else
        if abs(t - t_treat) < delta_t / 2
            p_1 = p_1 - mi * p_1 * delta_t + p1;
            m = 1;
        else
            p_1 = p_1 - mi * p_1 * delta_t;
        end
    end
    
    if C == 0
        f = 0;
    else
        f = p * (1 + p_1) * (T / C) ^ q / (s + (T / C) ^ q);
        if isnan(f)
            f = 0; %p * (1 + p_1);
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
%      Antigens natural release      %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = natural_release(rho, C)
    f = rho * C;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              %
%      T-cells activation      %
%                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [T, A, T_ac, T_in, m, c_4] = A_activate_T(a, b, K, h, c4, c_4, ni, t_treat, t, delta_t, T, A, vol_flag, time_flag)   
    m = 0;
    if vol_flag == 0 || time_flag == 0
    else
        if abs(t - t_treat) < delta_t / 2
            c_4 = c_4 - ni * c_4 * delta_t + c4;
            m = 1;
        else
            c_4 = c_4 - ni * c_4 * delta_t;
        end
    end

    T_ac = a * T * A;            % active
    T_in = b /(1 + c_4) * T * A; % inactive
    
    T0_flag = T + delta_t * (- T_ac - T_in + h * K) < 0; % T(t+1) < 0
    A0_flag = A + delta_t * (- T_ac - T_in) < 0; % A(t+1) < 0
    
    if T0_flag || A0_flag        
        if T0_flag
            delta_t_1 = -T / (- T_ac - T_in + h * K); %T = 0
            T_1 = 0;
            A_1 = max(0, A + delta_t_1 * (- T_ac - T_in));
            T_ac = delta_t_1 * T_ac;
            T_in = delta_t_1 * T_in;
        elseif A0_flag
            delta_t_1 = -A / (- T_ac - T_in); %A = 0
            A_1 = 0;
            T_1 = max(0, T + delta_t_1 * (- T_ac - T_in + h * K));
            T_ac = delta_t_1 * T_ac;
            T_in = delta_t_1 * T_in;
        else
            delta_t_2 = -A / (- T_ac - T_in); %A = 0
            delta_t_3 = -T / (- T_ac - T_in + h * K); %T = 0
            delta_t_1 = min(delta_t_2, delta_t_3);
            A_1 = 0;
            T_1 = 0;
            T_ac = delta_t_1 * T_ac;
            T_in = delta_t_1 * T_in;
            delta_t_1 = delta_t_3;
        end       
        delta_t_2 = delta_t - delta_t_1;
        T = max(0, min(K, T_1 + delta_t_2 * h * K));
        A = A_1;
    else
        T = max(0, min(K, T + delta_t * (- T_ac - T_in + h * K)));
        A = max(0, A + delta_t * (- T_ac - T_in));
        T_ac = delta_t * T_ac;
        T_in = delta_t * T_in;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                              %
%      Antigens release produced after RT      %
%                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = RT_release(psi, C)
    f = max(0, psi * C);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                %
%      Antigens leaving the system naturally     %
%                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = A_natural_out(sigma, A)
    f = - sigma * A;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                               %
%      T-cells leaving the system naturally     %
%                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = T_natural_out(eta, T)
    f = - eta * T;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                 %
%      T-cells leaving the system because of the interaction with tumor cells     %
%                                                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = immune_death_T(iota, C, T)
    f = - iota * T * C;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                    %
%      Tumor cell death kinetic      %
%                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = tum_kinetic(phi, tau_1, tau_2, t)
    if t <= tau_1
        a = 0;
    elseif t > tau_2
            a = 1;
    else
            a = (t - tau_1) / (tau_2 - tau_1);
    end

    f = - a * phi;
end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              %
%      Tumor volume (mm^3)     %
%                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = tumor_volume(C, T, vol_C, vol_T)
    f = C * vol_C + T * vol_T;
end


%%
%%%%%%%%%%%%%%%%%%%
%                 %
%      Markov     %
%                 %
%%%%%%%%%%%%%%%%%%%

function C = markov_TCP_analysis(im_death, prol, C, delta_t)
    cell_num = round(C);

    f = prol * delta_t;     %Birth probability
    g = im_death * delta_t; %Dead probability

    e = f + g;

    if e > 1
        f = f/e;
        g = g/e;
    end
    
    cell_array = randsrc(1, cell_num, [2, 0, 1; f, g, min( max(0, 1 - f - g), 1)]);
    C = sum(cell_array);
end